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Abstract 

The problem of error growth due to the incomplete knowledge of the evolution 
law which rules the dynamics of a given physical system is addressed. Major 
interest is devoted to the analysis of error amplification in systems with many 
characteristic times and scales. The importance of a proper parameterization 
of fast scales in systems with many strongly interacting degrees of freedom is 
highlighted and its consequences for the modelization of geophysical systems 
are discussed. 

I. INTRODUCTION 

The ability to predict the future state of a system, given its present state, stands at 
the foundations of scientific knowledge with relevant implications from an applicative point 
of view in geophysical and astronomical sciences. In the prediction of the evolution of a 
system, e.g. the atmosphere, we are severely limited by the fact that we do not know 
with arbitrary accuracy the evolution equations and the initial conditions of the system. 
Indeed, one integrates a mathematical model given by a finite number of equations. The 
initial condition, a point in the phase space of the model, is determined only with a finite 
resolution (i.e. by a finite number of observations) (Monin 1973). 

Using concepts of dynamical systems theory, there have been some progresses in under- 
standing the growth of an uncertainty during the time evolution. An infinitesimal initial 
uncertainty {5q — > 0) in the limit of long times {t —>■ oo) grows exponentially in time with 
a typical rate given by the leading Lyapunov exponent A, |5a;(t)| ~ 5oexp(At). Therefore if 
our purpose is to forecast the system within a tolerance A, the future state of the system 
can be predicted only up to the predictability time, given by: 
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In literature, the problem of predictability with respect to uncertainty on the initial condi- 
tions is referred to as predictability of the first kind. 

In addition, in real systems we must also cope with the lack of knowledge of the evolution 
equations. Let us consider a system described by a differential equation: 

^x(t) = f(x,t), x,fG7^^ (2) 

As a matter of fact we do not know exactly the equations, and we have to devise a model 
which is different from the true dynamics: 

^ x(t) = f,(x, t) where f,(x, t) = f (x, t) + e5f (x, t) . (3) 

Therefore, it is natural to wonder about the relation between the true evolution {reference 
or true trajectory xr(i:)) given by and that one effectively computed {perturbed or model 
trajectory XM{t)) given by (|^). This problem is referred to as predictability of the second 
kind. 

Let us make some general remarks. At the foundation of the second kind predictability 
problem there is the issue of structural stability (Guckenheimer et al. 1983): since the 
evolution laws are known only with finite precision it is highly desirable that at least certain 
properties are not too sensitive to the details of the equations of motion. For example, in 
a system with a strange attractor, small generic changes in the evolution laws should not 
change drastically the dynamics (see Appendix ^ for a simple example with non generic 
perturbation). 

In chaotic systems the effects of a small generic uncertainty on the evolution law are 
similar to those due to the finite precision on the initial condition (Crisanti et al. 1989). 
The model trajectory of the perturbed dynamics diverges exponentially from the reference 
one with a mean rate given by the Lyapunov exponent of the original system. The statistical 
properties (such as correlation functions and temporal averages) are not strongly modified. 
This last feature has been frequently related to the shadowing lemma (Guckenheimer et 
al. 1983; Ott 1993): almost all trajectories of the true system can be approximated by 
a trajectory of the perturbed system starting from a slightly different initial condition. 
However, as far as we know, the shadowing lemma can be proven only in special cases and 
therefore it cannot be straightforwardly invoked to explain the statistical reproducibility 
in a generic case. In addition, in real systems the size of an uncertainty on the evolution 
equations is determinable only a posteriori, based on the ability of the model equations to 
reproduce some of the features of the phenomenon. 

In dynamical systems theory, the problems of first and second kind predictability is 
essentially understood in the limit of infinitesimal perturbations. However even in this limit 
we must also consider the fluctuations of the rate of expansion which can lead to relevant 
modifications of the predictability time (0), in particular for strongly intermittent systems 
(Benzi et al. 1985; Paladin et al. 1987; Crisanti et al. 1993). 

As far as finite perturbations are considered, the leading Lyapunov exponent is not rele- 
vant for the predictability issue. In presence of many characteristic times and spatial scales 
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the Lyapunov exponent is related to the growth of small scale perturbations which saturates 
on short times and has very little relevance for the growth of large scale perturbations (Leith 
and Kraichnan 1972; Monin 1973; Lorenz 1996). To overcome this shortcoming, a suitable 
characterization of the growth of non infinitesimal perturbations, in terms of the Finite Size 
Lyapunov Exponent (FSLE), has been recently introduced (Aurell et al. 1996 and 1997). 

Also in the case of second kind predictability one has often to deal with errors which are 
far from being infinitesimal. Typical examples are systems described by partial differential 
equations (e.g. turbulence, atmospheric flows). The study of these systems is performed 
by using a numerical model with unavoidable severe approximations, the most relevant of 
which is the necessity to cut some degrees of freedom off; basically, the small scale variables. 

The aim of this paper is to analyze the effects of limited resolution on the large scale 
features. This raises two problems: in first place one has to deal with perturbations of 
the evolution equations which in general cannot be considered small; second, the parame- 
terization of the unresolved modes can be a subtle point. We shall show that the Finite 
Size Lyapunov Exponent is able to characterize the effects of uncertainty on the evolution 
laws. Moreover we shall discuss the typical difficulties arising in the parameterization of the 
unresolved scales. 

This paper is organized as follows. In section ^ we report some known results about the 
predictability problem of the second kind and recall the definition of the FSLE. In section 
m we present numerical results on simple models. In section |V| we consider more complex 
systems with many characteristic times. Section ^ is devoted to summarize the results. In 
Appendix ^ we illustrate a simple example of structural unstable system. In Appendix 
we describe the method for the computation of the FSLE and in Appendix y we discuss the 
problem of the parameterization of the unresolved variables. 



II. EFFECTS OF A SMALL UNCERTAINTY ON THE EVOLUTION LAW 

In the second kind predictability problem, we can distinguish three general cases depend- 
ing on the original dynamics. In particular, equation @ may display: 

(i) trivial attractors: asymptotically stable fixed points or attracting periodic orbits; 

(ii) marginally stable fixed points or periodic/quasi-periodic orbits as in integrable Hamil- 
tonian systems; 

(iii) chaotic behavior. 

In case (i) small changes in the equations of motion do not modify the qualitative features 
of the dynamics. Case (ii) is not generic and the outcome strongly depends on the specific 
perturbation 6i, i.e. it is not structurally stable. In the chaotic case (iii) one expects that 
the perturbed dynamics is still chaotic. In this paper we will consider only this latter case. 

Let us also mention that, in numerical computations of evolution equations (e.g. differen- 
tial equations), there are two unavoidable sources of errors: the finite precision representation 
of the numbers which causes the computer phase space to be necessarily discrete and the 
round-off which introduces a sort of noise. Because of the discrete nature of the phase space 
of the system studied on computer, orbits numerically computed have to be periodic. Nev- 
ertheless the period is usually very large, apart for very low computer precision (Crisanti et. 
al 1989). We do not consider here this source of difficulties. The round-off produces on eq. 
(0) a perturbation which can be written as 5f (x, t) = w(x) f (x, t) and e ~ 10~" [a =number 
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of digits in floating point representation) where w = 0(1) is an unknown function which 
may depend on f and on the software of the computer (Knut 1969). In general, the round-off 
error is very small and may have, as much as the noise, a positive role, as underlined by 
Ruelle (1979), in selecting the physical probability measure, the so-called natural measure, 
from the set of ergodic invariant measures. 

In chaotic systems the effects of a small uncertainty on the evolution law is, for many 
aspects, similar to those due to imperfect knowledge of initial conditions. This can be 
understood by the following example. Consider the Lorenz equations (Lorenz 1963) 

a{y - x) 

Rx ~ y — xz (4) 
xy — hz. 



dx 

Itt 
dy 

dt 

dz 



In order to mimic an experimental error in the determination of the evolution law we consider 
a small error e on the parameter R: R ^ R + e. Let us consider the difference Ax(t) = 
XA/(i) — ^xit) with, for simplicity, Ax(0) = 0, i.e. we assume a perfect knowledge of the 
initial conditions. One has, with obvious notation: 

dAx (9f df, 

^ = f.(xM)-f(x.)^-Ax+-e. (5) 

At time t = one has |Ax(0)| = 0, therefore |Ax(t)| grows initially only by the effect of the 
second term in (^. At later times, when |Ax(i)| ^ 0(e) the leading term of (^ becomes 
the flrst one, and we recover the flrst kind predictability problem for an initial uncertainty 
6o ~ e. Therefore, apart from an initial (not particularly interesting) growth, which depends 
strongly on the speciflc perturbation, the evolution of < log(|Ax(t)|) > follows the usual 
linear growth with the slope given by the leading Lyapunov exponent. Typically the value 
of the Lyapunov exponent computed by using the model dynamics differs from the true one 
by a small amount of order e, i.e. \m = Xt + 0(e) (Crisanti et al. 1989). 

This consideration applies only to inflnitesimal perturbations. The generalization to flnite 
perturbations requires the extension of the Lyapunov exponent to flnite errors. Let us now 
introduce the Finite Size Lyapunov Exponent for the predictability of flnite perturbations. 
The deflnition of FSLE X{6) is given in terms of the "doubling time" Tr{6), that is the time 
a perturbation of initial size 6 takes to grow by a factor r (> 1): 

where (■ ■ ■)( denotes average with respect to the natural measure, i.e. along the trajectory 
(see Appendix 0). For chaotic systems, in the limit of inflnitesimal perturbations {6 — > 0) 
X{6) is nothing but the leading Lyapunov exponent A (Benettin et al. 1980). Let us note 
that the above deflnition of X{6) is not appropriate to discriminate cases with A = and 
A < 0, since the predictability time is positive by deflnition. Nevertheless this is not a 
limitation as long as we deal with chaotic systems. 
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In many realistic situations the error growth for infinitesimal perturbations is dominated 
by the fastest scales, which are typically the smallest ones (e.g. small scale turbulence). 
When 6 is no longer infinitesimal, X{6) is given by the fully nonlinear evolution of the 
perturbation. In general A (5) < A, according to the intuitive picture that large scales 
are more predictable. Outside the range of scales in which the error 6 can be considered 
infinitesimal, the function A (5) depends on the details of the dynamics and in principle on 
the norm used. In fully developed turbulence one has the universal law A (5) ~ in the 
inertial range (Aurell et al. 1996 and 1997). It is remarkable that this prediction, which can 
be obtained within the multifractal model for turbulence, is not affected by intermittency 
and it gives the law originally proposed by Lorenz (1969). The behavior of X{6) as a function 
of 6 gives important information on the characteristic times and scales of the system and it 
has been also applied to passive transport in closed basins (Artale et al. 1997). 

Let us now return to the example (H). We compute \tt{S), the FSLE for the true 
equations, and Atm(^), the FSLE computed following the distance between one true tra- 
jectory and one model trajectory starting at the same point. These are shown in Figure 
0. The true FSLE Att(5) displays a plateau indicating a chaotic dynamics with leading 
Lyapunov exponent A ^ 1. Concerning the second kind predictability, for 6 > e the second 
term in (^ becomes negligible and we observe the transition to the Lyapunov exponent 
Xtm{S) — Att(^) — A. In this range of errors the model system recovers the intrinsic pre- 
dictability of the true system. For very small errors, 6 < e, where the growth of the error is 
dominated by the second term in (|^), we have Atm('5) > Att('^)- 

This example shows that it is possible to recover the intrinsic predictability of a chaotic 
system even in presence of some uncertainty in the model equations. 

The relevance of the above example is however limited by the fact that @ does not 
involve different scales. To investigate the effect of spatial resolution on predictability let 
us consider the advection of Lagrangian tracers in a given Eulerian field. We study a time- 
dependent, two dimensional, velocity field given by the superposition of large scale (resolved) 
eddies and small scale (possibly unresolved) eddies. 

The streamfunction we consider is a slight modification of a model originally proposed 
for chaotic advection in Rayleigh-Benard convection (Solomon et al. 1988): 

"^{x, y, t) = ^{x, y, t; ki, ujl, Bl) + e ■ ^{x, y, t; ks, ujs, Bs) (7) 

with 

ip{x, y, t; k, u,B) = — sin {k[x + B sm{ut)]} sin (ky) (8) 
k 

The first term represents the large-scale fiow, i.e. the resolved part of the fiow, the 
second one mimics the unresolved small scale term and e measures the relative amplitude. 
We choose ks ^ k^ and ojs ^ ojl m. order to have a sharp separation of space and time 
scales. 

The Lagrangian tracers evolve according to the equations: 

dx difj dy dip , . 

dt dy ' dt dx ' 

We use the complete stream function (0) for the true dynamics and only the large-scale term 
for the model dynamics. 
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The time dependence induces chaotic motion and diffusion in the x direction, without 
inserting any noise term (Solomon et al. 1988). For what concerns \'j't 

one observes three 

regimes (Figure ^). For very small errors 5 < 271 /ks the exponential separation is ruled by 
the fastest scale and \tt — A, i.e. we recover the Lyapunov exponent of the system. At 
intermediate errors we observe a second small plateau corresponding to the large-scale term 
for the model dynamics. For larger errors, 6 > 271 /ki one has \tt ~ i-e. diffusive 
behavior (see Artale et. al 1997). 

The model FSLE \tm{^) cannot recover the small scale features: for 5 <^ e we observe 
the scaling \tm{^) ~ which can be understood by the following argument. In this region 
the distance between the reference and the true trajectories grows as d5/dt ~ e and thus, 
by a dimensional estimate, one has: 

Atm(5) - ^ ~ ^ ■ (10) 

Nevertheless, for larger 5 the model fairly captures the small plateau displayed by \tt{^) 
which corresponds to the slow time scale, then at large 5 (i.e. for 5 greater than the 5 > 
2vr/fci) we recover the diffusive behavior with the correct diffusion coefficient. This last 
feature can be understood by the fact that the diffusion coefficient, being an asymptotic 
quantity of the flow, is not influenced by the details of small scale structures. 

This example is rather simple: large scales do not interact with the small ones and the 
number of degrees of freedom is very small. Therefore in this case the crude elimination of 
the small scale component does not prevent the possibility of a fair description of large scale 
features. 

In the following we will consider more complex situations, in which strongly interacting 
degrees of freedom with different characteristic times are involved. In these cases the correct 
parameterization of the unresolved modes is crucial for the prediction of large scale behavior. 



III. SYSTEMS WITH TWO TIME SCALES 

Before analyzing in detail the effects of non infinitesimal perturbations of the evolution 
laws in some specific models let us clarify our aims. We consider a dynamical system written 
in the following form: 

-rr = f(x,y) 

dt .^^N 

dy 

where f , x G 7^" and g, y G 71"^, in general n ^ m. Now, let us suppose that the fast variables 
y cannot be resolved: a typical example are the subgrid modes in PDE discretizations. In 
this framework, a natural question is: how must we parameterize the unresolved modes (y) 
in order to predict the resolved modes (x) ? 

As discussed by Lorenz (1996), to reproduce - at a qualitative level - a given phenomenol- 
ogy, e.g. the ENSO phenomenon, one can drop out the small scale features without negative 
consequences. But one unavoidably fails in forecasting the ENSO (i.e. the actual trajectory) 
without taking into account in a suitable way the small scale contributions. 
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An example in which it is relatively simple to develop a model for the fast modes is 
represented by skew systems: 

-77 = f(x,y) 

dt /^2) 

Tt = 

In this case, the fast modes (y) do not depend on the slow ones (x). One can expect that 
in this case, neglecting the fast variables or parameterizing them with a suitable stochastic 
process, should not drastically affect the prediction of the slow variables (Boffetta et al. 
1996). 

On the other hand, if y feels some feedback from x, we cannot simply neglect the 
unresolved modes. In Appendix ^ we discuss this point in detail. In practice one has 
to construct an effective equation for the resolved variables: 

^ = fM(x,y(x)), (13) 

where the functional form of y(x) and fM are found by phenomenological arguments and/or 
by numerical studies of the full dynamics. 

Let us now investigate an example with a recently introduced toy model of the atmo- 
sphere circulation (Lorenz 1996; Lorenz et al. 1998) including large scales Xk (synoptic 
scales) and small scales yj^k (convective scales): 



dxk 
dt 

dyj,k 
dt 



k 



-Xk-l {Xk-2 - Xk+l) -UXk + F - Y.j=l Vj 

(14) 

-cbyj+i^k {yj+2,k - yj~i,k) - cuyj^k + Xk 



where k = 1,...,K and j = 1,...,J. As in (Lorenz 1996) we assume periodic boundary 
conditions on k {xx+k = Xk, yj,K+k = yj,k) while for j we impose yj+j,k = yj,k+i- The 
variables Xk represent some large scale atmospheric quantities in K sectors extending on a 
latitude circle, while the yj^k represent quantities on smaller scales in J ■ i^' sectors. The 
parameter c is the ratio between fast and slow characteristic times and b measures the 
relative amplitude. 

As pointed out by Lorenz, this model shares some basic properties with more realistic 
models of the atmosphere. In particular, the non-linear terms, which model the advection, 
are quadratic and conserve the total kinetic energy J^ki^l + J^j y'j^k) ^^e unforced (F = 0), 
inviscid {u = 0) limit; the linear terms containing u mimic dissipation and the constant term 
F acts as an external forcing preventing the total energy from decaying. 

If one is interested in forecasting the large scale behavior of the atmosphere by using 
only the slow variables, a natural choice for the model equations is: 

dxk 

— = -Xk-l ixk-2 - Xk+i) -uxk + F - G'fc(x) , (15) 

where (^^(x) represents the parameterization of the fast components in (|T4D (see Appendix 
0). 
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The FSLE for the true system (Boffetta et al. 1998) is shown in Figure ^ and displays 
the two characteristic plateau corresponding to fast component (for 6 -C 0.1) and slow com- 
ponent for large 6 . Figure |^ also shows what happens when one simply neglects the fast 
components yj^k (i-e. G(x) = 0). At very small 6 one has \tm{S) — as previously 
discussed. For large errors we observe that, with this rough approximation, we are not able 
to capture the characteristic predictability of the original system. More refined parameter- 
izations in terms of stochastic processes with the correct probability distribution function 
and correlation times do not improve the forecasting ability. 

The reason for this failure is due to the presence of a feedback term in the equations 
(p!^) which induces strong correlations between the variable and the unresolved coupling 
J2j=i yj,k- For a proper parameterization of the unresolved variables we follow the strategy 
discussed in Appendix |C[ Basically we adopt 

G(x) = UeXk , (16) 

in which is a numerically determined parameter. Figure ^ shows that, although small 
scale are not resolved, the large scale predictability is well reproduced and one has Xtm{S) — 
Xtt{S) for large 6. We conclude this section by observing that the proposed parameterization 
(p!6|) is a sort of eddy viscosity parameterization. 



IV. LARGE SCALE PREDICTABILITY IN A TURBULENCE MODEL 

We now consider a more complex system which mimics the energy cascade in fully 
developed turbulence. The model is in the class of the so called shell models introduced 
some years ago for a dynamical description of small-scale turbulence. For a recent review 
on shell models see Bohr et al. 1998. This model has relatively few degrees of freedom 
but involves many characteristic scales and times. The velocity field is assumed isotropic 
and it is decomposed on a finite set of complex velocity components u„ representing the 
typical turbulent velocity fluctuation on a "shell" of scales in = In order to reach 

very high Reynolds number with a moderate number of degrees of freedom, the scales are 
geometrically spaced as kn = ko2"' {n = 1, ...A^). 

The specific model here considered has the form (L'vov et al. 1998) 

^ = i " ^^n^^i^-lMn+l + ^fcn-lMn-2«n-l^ - l^klUn + fn (17) 

where u represent the kinematic viscosity and is a forcing term which is restricted only 
to the first two shells (in order to mimic large scale energy injection). 

Without entering in the details, we recall that the Shell Model (|T7D displays an energy 
cascade d la Kolmogorov from large scales (small n) to dissipative scales [n ~ A^) with 
a statistical stationary energy flux. Scaling laws for the average velocity components are 
observed: 

(Kl) ^ K^' (18) 

with exponents close to the Kolmogorov 1941 values Cp = p/3. 
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From a dynamical point of view, model ([T7|) displays complex chaotic behavior which is 
responsible of the small deviation of the scaling exponents (intermittency) with respect to 
the Kolmogorov values. Neglecting this (small) intermittency effects, a dimensional estimate 
of the characteristic time (eddy turnover time) for scale n gives 



Ur. 



K'" ■ (19) 



The scaling behavior holds up to the Kolmogorov scale 7] = l/k^ defined as the scale 
at which the dissipative term in (|l^ becomes relevant. The Lyapunov exponent of the 
turbulence model can be estimated as the fastest characteristic time and one has the 
prediction (Ruelle 1979) 

A ^ 1 ^ Re^/^ (20) 

where we have introduced the Reynolds number Re oc 1/z/. It is possible to predict the 
behavior of the FSLE by observing that the faster scale A;„ at which an error of size S is 
still active (i.e. below the saturation) is such that m„ ~ 6. Thus A(5) ~ l/r„ and, using 
Kolmogorov scaling, one obtains (Aurell et al. 1996 and 1997) 

Att(^) - I U T ^ -<1 < (21) 
^ ' [ 5 for < () < Mo 

To be more precise there is an intermediate range between the two showed in (|21|). For a 
discussion on this point see (Aurell et al. 1996 and 1997). 

In order to simulate a finite resolution in the model, we consider a modelization of ( P!7[ ) 
in terms of an eddy viscosity (Benzi et al. 1998) 



dt 



kn+lU*^+^Un+2 - ^/c„<_iUn+l + ]^kn-lUn-2Un-l^ " Z^i^^/c^Mn + fn (22) 



where now n = 1, ...,Nm < N and the eddy viscosity, restricted to the last two shells, has 
the form 



^ = f^-T^ «N,.i-l + 6n,NM) (23) 



where k is a constant of order 1 (see Appendix |CD. The model equations (^2]) are the 
analogous of large eddy simulation (LES) in Shell Model which is one of the most popular 
numerical method for integrating large scale flows. Thus, although Shell Models are not 
realistic models for large scale geophysical flows (being nevertheless a good model for small 
scale turbulent fluctuations), the study of the effect of truncation in term of eddy viscosity 
is of general interest. 

In Figure |^ we show \mm{S), i-e. the FSLE computed for the model equations ( P^] ) with 
= 24 at different resolutions Nm = 9, 15, 20. A plateau is detected for small amplitudes of 
the error 6, corresponding to the leading Lyapunov exponent, which increases with increasing 
resolution - being proportional to the fastest timescale - according to A ~ ^n^' larger 
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S the curves collapse onto the Xtt{S), showing that large-scale statistics of the model is not 
affected by the small-scales resolution. 

The capability of the model to predict satisfactorily the statistical features of the "true" 
dynamics is not anyway determined by XMMi^) but by Atm(^), which is shown in Figure 

Increasing the resolution Nm = 9, 15, 20 towards the fully resolved case N = 24 the 
model improves, in agreement with the expectation that Atm approaches Xtt for a perfect 
model. At large 6 the curves practically coincide, showing that the predictability time for 
large error sizes (associated with large scales) is independent on the details of small-scale 
modeling. Better resolved models achieve Atm — Xtt for smaller values of the error 6. 

V. CONCLUSIONS 

In this Paper the effects of the uncertainty of the evolution laws on the predictability 
properties are investigated and quantitatively characterized by means of the Finite Size Lya- 
punov Exponent. In particular, we have considered systems involving several characteristic 
scales and times. In these cases, it is rather natural to investigate what is the effect of small 
scale parameterization on large scale dynamics. 

It has been shown that in systems where there is a negligible feedback on the small scales 
by the large ones, the dynamics of the former ones can be thoroughly discarded, without 
affecting the statistical features of large scales and the ability to forecast them. On the 
other side, when this feedback is present, the crude approximation of cutting the small scale 
variables off is no longer acceptable. In this case one has to model the action of fast modes 
(small scales) on slow modes (large scales) with some effective term, in order to recover 
a satisfactory forecasting of large scales. The renowned eddy-viscosity modelization is an 
instance of the general modeling scheme that has been here discussed. 
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APPENDIX A: AN EXAMPLE OF STRUCTURAL UNSTABLE SYSTEM 

In order to see that a non generic perturbation, although very "small", can produce 
dramatic changes in the dynamics, let us discuss a simple example following (Berkooz 1994; 
Holmes et al. 1996). We consider the one-dimensional chaotic map Xt+i = f{xt) with 
f{x) = 4x mod 1, and a perturbed version of it: 
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2^+3 



X G 



247 265 
384' 384 



(Al) 



6 ^ 

4x mod 1 otherwise . 



265 17 
384' 24 



The perturbed map is identical to the original outside the interval [5/8,17/24], and the 
perturbation is very small in L2 norm. Nevertheless, the fixed point a; = |, which is 
unstable in the original dynamics, becomes stable in the perturbed one. Moreover it is 
a global attractor for fp{x), i.e. almost every point in [0, 1] asymptotically approaches a; = | 
(see Figure |l|). 

Now, if one compares the trajectories obtained iterating f[x) or fp{x) it is not difficult 
to understand that orbits starting outside [5/8, 17/24] remain identical for a certain time 
but unavoidably they differ utterly in the long time behavior. It is easy to realize that 
the transient chaotic behavior of the perturbed orbits can be rendered arbitrarily long by 
reducing the interval in which the two dynamics differ. This example shows how even an 
ostensibly small perturbation (in usual norms) can modify dramatically the dynamics. 



APPENDIX B: COMPUTATION OF THE FINITE SIZE LYAPUNOV 

EXPONENT 

In this appendix we discuss in detail the computation of the Finite Size Lyapunov Expo- 
nent for both continuous dynamics (differential equations) and discrete dynamics (maps). 

The practical method for computing the FSLE goes as follows. Defined a given norm for 
the distance 6{t) between the reference and perturbed trajectories, one has to define a series 
of thresholds 5„ = r"5o {n = 1, . . . , N), and to measure the "doubling times" Tr{6n) that a 
perturbation of size 6n takes to grow up to Sn+i- The threshold rate r should not be taken 
too large, because otherwise the error has to grow through different scales before reaching 
the next threshold. On the other hand, r cannot be too close to one, because otherwise the 
doubling time would be of the order of the time step in the integration. In our examples we 
typically use r = 2 or r = a/2. For simplicity is called "doubling time" even if r 7^ 2. 

The doubling times Tr{6n) are obtained by following the evolution of the separation from 
its initial size 6min ^ up to the largest threshold 6n- This is done by integrating the 
two trajectories of the system starting at an initial distance 5mm- In general, one must 
choose Smin So, in order to allow the direction of the initial perturbation to align with 
the most unstable direction in the phase-space. Moreover, one must pay attention to keep 

< ^saturation, SO that all the thrcsholds can be attained {Ssaturation is the typical distance of 
two uncorrelated trajectory, i.e. the size of the attractor). For the second kind predictability 
problem, i.e. the computation of \tm{S), one can safely take 6min = because this do not 
prevent the separation of trajectories. 

The evolution of the error from the initial value 6min to the largest threshold 6n carries 
out a single error-doubling experiment. At this point one rescales the model trajectory at 
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the initial distance 5mm with respect to the true trajectory and starts another experiment. 
After A/" error-doubhng experiments, we can estimate the expectation value of some quantity 
A as: 

{A)e = YrT.A. (Bl) 

■'^ i=i 

This is not the same as taking the time average as in (j^) because different error doubling 
experiments may takes different times. Indeed we have 

{Ah = U' Am = = iil^ . (B2) 

T Jo Ei Ti (r)e 



In the particular case in which A is the doubling time itself we have from and ( |B2|) 



A(5„) = --^lnr. (B3) 

{Ir[dn)}e 

The method described above assumes that the distance between the two trajectories is 
continuous in time. This is not true for maps of for discrete sampling in time and the method 
has to be slightly modified. In this case Tj.(5„) is defined as the minimum time at which 
5{Tr) > r6n. Because now S(Tr) is a fluctuating quantity, from (p^ ) we have 

We conclude by observing that the computation of the FSLE is not more expensive 
than the computation of the Lyapunov exponent by standard algorithm. One has simply to 
integrate two copies of the system (or two different systems for second kind predictability) 
and this can be done also for very complex simulations. 



APPENDIX C: PARAMETERIZATION OF SMALL SCALES 

Typically a realistic problem (e.g. turbulence) involves many interacting degrees of 
freedom with different characteristic times. Let us indicate with z the state of the system 
under consideration, with an evolution law: 

^ = F(z), F,zG7^^. (CI) 

The dynamical variables z can be split in two sets: 

z=(x,y), (C2) 

where x G TZ"' and y G TZ"^ {N = n + m), being respectively x and y the "slow" and "fast" 
variables. The distinction between slow and fast variables is often largely arbitrary. 

The evolution equation ( PT| ) is divided into two blocks, the first one containing the 
dynamics of the slow variables, the second one associated with the dynamics of the fast 
variables 
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^ = Fi(x) + F2(x,y) 
< ^* (C3) 
^^ = Fi(x,y) + F,(y) 

If one is interested only in the slow variables it is necessary to write an "effective" 
equation for x. As far as we know there is only one case for which it is simple to find the 
effective equations for x. If the characteristic times of the fast variables are much smaller 
than those ones of the x (adiabatic limit), one can write: 

y =< y > +77(t) (C4) 

where 77 is a Wiener process, i.e. a zero mean Gaussian process with 

< v^{t)r]j{t) >=< 6y^ > 6,,6{t - t) . (C5) 

Therefore one obtains for the slow variables: 

^ = Fi(x) + 5Fi(x) + <5W(x, 77) (C6) 

where 5Fi(x) = F2(x,<y>) + 5F2, SF^j = l/2Y.id^F^J dy^^ < 5yf > and SW^ = 
J2i 9F2J /dyjl^y^ Vi{t)- Basically the slow variables x obey to a non linear Langevin equa- 
tion. 

Here the role of the fast degrees of freedom becomes relatively simple: they give small 
changes to the drift Fi Fi + 5Fi and a noise term 5W(x, 77). We remark that the validity 
of the above argument is rather limited. Even if one has a large time scale separation, the 
statistics of the fast variables can be very far from the Gaussian distribution. In particular, 
in system with feedback (Fi 7^ 0) one cannot model the fast variable y independently of the 
resolved x. 

In the generic situation the construction of the effective equation for x requires to follow 
phenomenological arguments which depend on the physical mechanism of the particular 
problem. For example, for the Lorenz '96 model discussed in sect. |T|, where F2,fc(x, y) = 
J2j=i,jyj,k, we use the following procedure for the parameterization of the fast variables and 
the building of the effective eq. for x. Instead of assuming (|C^ ) we mimic the fast variables 
in terms of the slow ones: 



y(t)=g(x(t))=<y|x(t)>+r7(t) 



(C7) 



where < |x > stands for the conditional average and ri{t) is a noise term. Inserting (C7) 
into the first of (|C3|) one obtains 



dx 

dt 



Fi(x) + F2(x, y) = Fi(x) + F2(x, <y|x>) + 5F2(x) 



(C8) 



where 



= E 



d'^F 



2A 



j,k ^yj 



dvjdyk 



iVjVk) 



(C9) 



y=<y\x> 
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In the Lorenz '96 model (0), because of the hnear couphng between the different scales, 
the terms 6F2 are absent and one has a close model for the large scale variables 



dx 

— = Fi(x) + F2(x,<y|x>) (CIO) 



The ansatz (|C7|) is well verified in the numerical simulations. We have computed the 
^tm{S) by using a best fit for F2 and we have obtained a good reproduction of the Arr('5) 
for large S. In the Lorenz '96 model (p!4D, where the coupling between slow and fast variables 
is practically linear, one has that -F2,fe(x, <y|x>) = J2j=i,j <yj,k\xk>— J^e^k- 

Now we will discuss the case of the Shell Model parameterization which pertains to the 
general issue of the subgrid-scale modelization. The literature on this field and the related 
problems (e.g. closure in fully developed turbulence) is enormous and we do not pretend to 
discuss here in details this field. Let us only recall the basic idea introduced over a century 
ago by Boussinesq, and later developed further by Taylor, Prandtl and Heisenberg - to cite 
some of the most famous ones- for fully developed turbulence (Frisch 1995). In a nutshell 
the idea is to mimic the energy flux from the large to the small scales (in our terms from 
slow to fast variables) by an effective dissipation: the effect of the small scales on the large 
ones can be modeled as an enhanced molecular viscosity. 

By simple dimensional arguments one can argue that the effects of small scales can be 
replaced by an effective viscosity at scales r, given by 

i/^*^) ~ r6v{r) (Cll) 

where 6v{r) is the velocity fluctuation on the scale r. 

The above argument for the Shell Model (|T^ gives (Benzi et al. 1998): 

^(^) = n^-^ (C12) 



k 



n 



where k ~ 0(1) is an empirical constant. From eq. ( |(J11D one could naively think to use 
dimensional argument d la Kolmogorov to set a constant eddy viscosity z/^*^^ ~ ^n^^^- ^^^^ 
way one forgets the dynamics and this can cause numerical blow up. More sophisticated 
arguments that do not include the dynamics lead to similar problems. 

Let us remark that the parameterization (|C12|) is not exactly identical to those obtained 



by closure approaches where the eddy viscosity is given in terms of averaged quantities. In 
our case this would mean to write {uD^^"^ instead of in ( |C12|) . 



After this discussion it is easy to recognize that the parameterization in terms of con- 
ditional averages introduced for the Lorenz '96 model is, a posteriori, an eddy viscosity 
model. 
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FIGURE CAPTIONS 



FIGURE [T|: The map fp of equation (|A1|) (solid line) and the original chaotic map / (dashed line). 

FIGURE]^: Finite Size Lyapunov Exponents Xtt{S) (+) and \tm{S) (x) versus S for the Lorenz 
model (^ with a = c = 10, b = 8/3, i? = 45 and e = 0.001. The dashed line represents 
the leading Lyapunov exponent for the unperturbed system (A ~ 1.2). The statistics 
is over 10^ realizations. 

FIGURE ||: \tt{S) (crosses, x) and \tm{S) (open squares, □) versus 6 for the Rayleigh-Benard 
model (0) with C = 0.5, fc^ = ul = 1, Bl = 0.3, ks = A, us = 4, Bs = 0.3 
and e = 0.125. The straight line indicates the slope. The statistics is over 10"^ 
realizations. 

FIGURE ^: Finite Size Lyapunov Exponents for the Lorenz '96 model Xtt{^) (solid line) and 
Xtm{S) versus 6 obtained by dropping the fast modes (+) and with eddy viscosity 
parameterization (x) as discussed in ( |T5D and ([l6|). The parameters are F = 10, 
K = 36 , J = 10, z/ = 1 and c = 6 = 10, implying that the typical y variable is 10 
times faster and smaller than the x variable. The value of the parameter z/g = 4 is 
chosen after a numerical integration of the complete equations as discussed in Appendix 
|C[ The statistics is over 10^ realizations. 

FIGURE 1^: The FSLE for the eddy- viscosity shell model (|2^) Xmm{S) at various resolutions Nm = 
9(+), 15(x), 20(*). For comparison it is drawn the FSLE Xtt{S) (continuous line). 
Here n = 0.4, ko = 0.05. 

FIGURE 1^: The FSLE between the eddy- viscosity shell model and the full shell model Xtm{S), at 
various resolutions Nm = 9(+), 15(x), 20(*). For comparison it is drawn the FSLE 
Xtt{S) (continuous line). The total number of shell for the complete model is = 24, 
with ko = 0.05, z/ = 10"^. 
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